#include "math.h"
#include <algorithm>
#include <string>
#include <sstream>
#include <fstream>
#include <iostream>
#include <vector>
#include <map>
#include <sys/time.h>
#include <unistd.h>
#include "Matrix.h"
#include <limits>

using namespace Numeric_lib;
using namespace std;


#ifndef NW_WEIGHTS
#define NW_WEIGHTS
/*
подробнее см. Biological Sequence analysis, Richard Durbin, p. 55
*/
#define REWARD      4   // вознаграждение за совпадение символов
#define REWARD_PUPY 2   // вознаграждение за совпадение символов (приведённых к purine-pyrimedine)
#define PENALTY     8   // штраф за несовпадение символов
#define OGAP       20   // штраф за открытие разрыва
#define EGAP       10   // штраф за продолжение разрыва
#endif

#ifndef USE_PUPY_WEIGHTS
#define USE_PUPY_WEIGHTS false
#endif

// возвращает текущее время в формате sec.msec
double cur_time() {
  struct timeval tv;
  gettimeofday(&tv, 0);
  return (double)tv.tv_sec + (double)tv.tv_usec*0.000001;
}

struct MICell {
    int m;
    int i;
};

void print_matrix_m(MICell* matrix, int l1, int l2) {
    for (int i1=0; i1<=l1; i1++) {
        for (int i2=0; i2<=l2; i2++) {
            cout << matrix[i1*(l2+1)+i2].m << "\t";
        }
        cout << "\n";
    }
}

void print_koord(int l1, int l2, int index) {
    int i = index % (l2+1);
    int j = index / (l2+1);
    cout << "KOORD: i = " << i << "; j = " << j << "\n";
}

void print_matrix_i(MICell* matrix, int l1, int l2) {
    for (int i1=0; i1<=l1; i1++) {
        for (int i2=0; i2<=l2; i2++) {
            cout << matrix[i1*(l2+1)+i2].i << "\t";
        }
        cout << "\n";
    }
}

int _nw_score(const string& seq1, const string& seq2, Matrix<MICell, 2>& matrix) {
    int i,j,s;
    int l1 = seq1.size();
    int l2 = seq2.size();
    MICell cell1, cell2;
    
    // initializing
    matrix[0][0].m = 0;
    matrix[0][0].i = 0;
    for (i=1; i<=l1; i++) {
        matrix[i][0].m = - OGAP - (i - 1) * EGAP;
        matrix[i][0].i = - OGAP - (i - 1) * EGAP;
    }
    for (i=1; i<=l2; i++) {
        matrix[0][i].m = - OGAP - (i - 1) * EGAP;
        matrix[0][i].i = - OGAP - (i - 1) * EGAP;
    }
    
    // calculate score
    for (i=1; i<=l1; i++) {
        for (j=1; j<=l2; j++) {
            if (seq1[i-1] == seq2[j-1]) {
                s = REWARD;
            }
            else {
                s = - PENALTY;
            }
            cell1 = matrix[i-1][j-1];
            matrix[i][j].m = s + max(cell1.m, cell1.i);
            cell1 = matrix[i-1][j];
            cell2 = matrix[i][j-1];
            matrix[i][j].i = max(max(cell1.m-OGAP, cell1.i-EGAP), max(cell2.m-OGAP, cell2.i-EGAP));
        }
    }
    return max(matrix[l1][l2].m, matrix[l1][l2].i);
}

int _nw_score_restricted(const string& seq1, const string& seq2, MICell* matrix) {
    int i,j,s,t1,t2;
    int l1 = seq1.size();
    int l2 = seq2.size();
    int MIN_INT = numeric_limits<int>::min();
    MICell cell;
    int k = max(abs(l1-l2), min(l1,l2)/10);

    //cout << seq1 << "\n";
    //cout << seq2 << "\n";

    int max_i1, max_i2;
    
    // initializing
    matrix[0].m = 0;
    matrix[0].i = 0;
    for (i=1; i<=k; i++) {
        matrix[i*(l2+1)].m = - OGAP - (i - 1) * EGAP;
        matrix[i*(l2+1)].i = - OGAP - (i - 1) * EGAP;
        matrix[i].m = - OGAP - (i - 1) * EGAP;
        matrix[i].i = - OGAP - (i - 1) * EGAP;
    }

    // calculate score
    int tmp_index;
    for (i=1; i<=l1; i++) {
        t1 = max(1, i-k);
        t2 = min(l2, i+k);
        for (j=t1; j<=t2; j++) {
            if (USE_PUPY_WEIGHTS) {
                switch (abs(seq1[i-1] - seq2[j-1])) {
                    case 0:  // equal
                        s = REWARD;
                        break;
                    case 6:  // a-g or g-a
                    case 17: // c-t or t-c
                        s = REWARD_PUPY;
                        break;
                    default: // mismatch
                        s = - PENALTY;
                        
                }
            }
            else {
                if (seq1[i-1] == seq2[j-1])
                    s = REWARD;
                else
                    s = - PENALTY;
            }
            tmp_index = i*(l2+1)+j;
            cell = matrix[tmp_index-l2-2];
            matrix[tmp_index].m = s + max(cell.m, cell.i);
            //matrix[tmp_index].m = -99;
            
            if (j<i+k) {
                cell = matrix[tmp_index-l2-1];
                max_i1 = max(cell.m-OGAP, cell.i-EGAP);
            } else {
                max_i1 = MIN_INT;
            }
            //if ((i==2) and (j==8)) {
            //    cout << "STRANGE CELL = " << matrix[tmp_index].i << ", max_i1=" << max_i1 << ", max_i2=" << max_i2 << "\n";
            //    cout << "STRANGE CELL: open=" << cell.m-OGAP << ", extend=" << cell.i-EGAP << "; MAX=" << max(cell.m-OGAP, cell.i-EGAP) <<"\n";
            //    cout << "STRANGE CELL: prev.m=" << cell.m << ", prev.i=" << cell.i << "\n";
            //    print_koord(l1, l2, tmp_index);
            //    print_koord(l1, l2, tmp_index-l2-2);
            //    print_koord(l1, l2, tmp_index-l2-1);
            //    print_koord(l1, l2, tmp_index-1);
            //}
            if (j>i-k) {
                cell = matrix[tmp_index-1];
                max_i2 = max(cell.m-OGAP, cell.i-EGAP);
            } else {
                max_i2 = MIN_INT;
            }
            matrix[tmp_index].i = max(max_i1, max_i2);
            //matrix[tmp_index].i = -99;
            
        }
    }
    //print_matrix_i(matrix, l1, l2);
    //printf("\n");
    //print_matrix_m(matrix, l1, l2);
    return max(matrix[l1*(l2+1)+l2].m, matrix[l1*(l2+1)+l2].i);
}

double nw_score(string seq1, string seq2) {
    int l1 = seq1.size();
    int l2 = seq2.size();
    Matrix<MICell, 2> matrix(l1+1, l2+1);
    
    return _nw_score(seq1, seq2, matrix);
}


double nw_score_mk(string seq1, string seq2,
                   MICell* matrix,
                   int iter_num) {
    /*
    seq1, seq2 - выравниваемые последовательности
    matrix - выделенная структура в памяти для матрицы замен
    iter_num - число итераций в Монте-Карло
    */
    
    int i;
    int real_score;
    int rand_score;
    double scores_sum = 0;
    vector<int> scores(iter_num);

    real_score = _nw_score_restricted(seq1, seq2, matrix);
    // Monte-Carlo
    for (i=0; i<iter_num; i++) {
        random_shuffle(seq1.begin(), seq1.end());
        rand_score = _nw_score_restricted(seq1, seq2, matrix);
        scores_sum += rand_score;
        scores[i] = rand_score;
        //cout << "MK iter=" << i << "\t score=" << rand_score << "\n";
    }
    
    // calculating Z
    double ev = scores_sum / iter_num;
    double diff_sum = 0;
    for (i=0; i<iter_num; i++) { diff_sum += (scores[i] - ev)*(scores[i] - ev); }
    double sd = sqrt(diff_sum/iter_num);

    //cout << "MK preview: " << real_score << "\t" << ev << "\t" << sd << "\t" << seq1.length() << "\n";
    //sort(scores.begin(), scores.end());
    //cout << "MK preview: " << scores[0] << "\t" << scores[iter_num/100*95] << "\t" << scores[iter_num-1] << "\n";

    if (sd > 0) {
        return (real_score - ev) / sd;
    } else {
        return 0.0;
    }
}
